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shocked  medium  bounded  by  the  contact  surface,  wedge/and  reflected  shock.  In  region 
(b)  die  solution  is  known  immediately  in  terms  of^Mj,  6^,  and  the  conditions  in  (a). 

Hie  shapes  of  the  Mach  stem,  contact  surface  and  bow  shock  are  expressed  parame 
as  I  -  F(s),  17  -  G(s).  Then  pc,  uc,  vc 
variables  in  double  power  series,  e.g., 


5U 

pc,  and  pd,  ud,  vd,  pd  are  obtained  by  expanding 


Pc(f,  17)  -  2  pfiV. 

id  <1 

substituting  in  the  ideal  fluid  equations,  and  equating  coefficients  of  like  powers  through 
some  order  N  *  max(i+j).-^The  resulting  algebraic  equations  are  solved  subject  to  the 
additional  relations  obtained  by  applying  the  reflection  conditions  on  the  wedge,  together 
with  tiie  jump  conditions  on  the  boundaries  ac  and  bd,  approximated  by  power  series 
expansions  of  the  F  and  G  functions.  Since  all  these  equations  sre  nonlinear,  solutions 
are  obtained  by  iteration  with  N  increasing  until  convergence  is  obtained.  The  Ben-Dor 
equation  for  the  fluid  quantities  in  regions  c,  d  at  the  triple  point  is  used  to  give  initial 
values.  -Becaiueiariation  within  each  region  is  smooth,  effectively  exact  descriptions  of 
most  features  of  interest  can  be  obtained  using  series  with  £  20  terms.  There  are  thus 
^  200  quantities  in  the  discretizatio&jifUfcc'problem,  compared  with  ^  104  in  a  con¬ 
ventional  finite-difference  treatment.^  The  method  generalizes  readily  to  complex  and 
double  Mach  reflections. 
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INTRODUCTION 

The  problem  of  a  planar  shock  reflecting  from  an  oblique  surface  goes 
back  over  a  hundred  years  to  Ernst  Mach.  Although  this  problem  is  important 
in  its  own  right,  much  of  the  interest  in  it  arises  because  of  the  need  for 
better  understanding  of  Mach  reflection  in  more  complicated  situations.  The 
field  has  been  the  object  of  particular  interest  during  the  last  thirty 
years;  the  experimental  and  theoretical  research  carried  out  during  this 
period  have  been  reviewed  by  Ben-Dor.1 

A  constant  planar  shock  propagating  into  a  uniform  ambient  gas  gives 
rise  in  the  absence  of  reflection  to  a  second  medium  with  uniform  thermody¬ 
namic  properties  in  the  region  behind  the  shock .  If  it  propagates  in  a  shock 
Manuscript  approved  September  13,  1983. 


tube  whose  walls  are  not  parallel  to  the  direction  of  propagation  (because  a 
wedge  has  been  inserted  along  the  side),  a  reflected  shock  wave  propagates 
back  into  the  interior  of  the  shock  tube.  When  the  wedge  angle  is  large,  so 
that  the  primary  shock  is  incident  nearly  normally  on  it,  the  reflected  shock 
is  also  planar  and  no  other  gasdynamic  discontinuities  appear  near  the 
reflection  point.  As  the  wedge  angle  decreases,  so  that  the  incident  shock 
becomes  more  and  more  nearly  glancing,  it  becomes  impossible  for  a  fluid 
particle  to  traverse  both  incident  and  reflected  shocks  and  still  "turn  the 
corner”  enough  to  end  up  moving  parallel  to  the  wedge  surface.  At  least  one 
additional  shock  (the  Mach  stem)  must  appear  (Fig.  1),  intersecting  the 
others  at  a  so-called  triple  point.  Because  some  of  the  material  in  the  zone 
between  the  Mach  stem  and  the  reflected  wave  has  been  shocked  once  and  some 
twice,  another  gasdynamic  discontinuity  (a  contact  surface)  must  also  extend 
from  the  triple  point  into  this  region,  terminating  somewhere  on  the  wedge 
surface.  Hie  reflected  shock  may  terminate  at  the  corner  of  the  wedge 
(attached  shock),  or  upstream  from  this  point  (detached  shock),  or  may  run 
into  a  second  triple  point  between  the  first  one  and  the  corner  (double  Mach 
reflection) .  The  latter  case  occurs  in  general  for  smaller  wedge  angles  than 
does  single  Mach  reflection;  an  intermediate  case  (complex  Mach  reflection) 
is  also  observed. 

One  would  like  to  derive  a  theoretical  description  of  Mach  reflection 
which  would  complement  the  experimental  results  and  address  some  of  the 
questions  the  latter  leave  unanswered,  such  as  the  structure  of  the  contact 
surface  near  the  wedge  surface,  and  whether  a  "triple  Mach”  regime  exists. 

An  analytical  solution  is  out  of  the  question,  though  pieces  of  the  problem 
(e.g.,  the  flow  in  the  neighborhood  of  the  triple  point1)  can  be  solved. 
Recourse  must  therefore  be  had  to  numerical  methods. 


Figure  1 « 

ew  -  io°. 


Interferogram  of  single  Mach  reflection  in  U2  (Mg  »  4.7 2 , 
P_  «  15  torr)  after  Ben-Dor1,  with  X  and  Y  axes  and 

Cl 


gasdynamic  discontinuities  drawn  in,  and  regions  labeled  with  a,  b,  c,  d 


Straightforward  differencing  of  the  fluid  equations  (using  either  an 
ideal  or  realistic  equation  of  state)  has  been  carried  out  with  considerable 
success1’'2.  A  rectilinear  Eulerian  mesh  with  M 00  -  200  zones  along  each 
axis,  possibly  varying  in  size  as  a  function  of  position  to  improve  resolu¬ 
tion  near  the  Mach  stem,  is  used.  The  conditions  ahead  of  and  behind  the 
incident  shock  are  used  as  boundary  conditions,  along  with  a  reflection  con¬ 
dition  on  the  wedge.  State-of-the  art  shock  capturing  techniques  such  as 
Flux-Corrected  Transport  ( FCT) 3  resolve  the  envelope  (reflected)  waveform 
accurately  and  permit  almost  all  the  other  gasdynamic  discontinuties  in  the 
problem  to  be  distinguished  (Fig.  2).  To  date,  however,  such  numerical 
solutions  have  not  surpassed  experimental  interferometric  data  in  accuracy, 
nor  have  they  succeeded  in  answering  any  of  the  outstanding  questions 
associated  with  Mach  reflection. 

Moreover,  such  calculations  leave  a  distinct  impression  of  brute  force. 
Advancing  the  four  fluid  equation  10 3  -  104  times teps  on  10 4  or  more  zones 
seems  profligate,  particularly  in  view  of  the  smoothness  of  most  of  the  gas- 
dynamic  discontinuities  and  the  gentle  variation  of  the  fluid  quantities  in 
the  regions  they  bound.  In  much  of  the  domain  of  the  calculation  the  solu¬ 
tion  is  known  a  priori  and  does  not  change  in  time. 

Furthermore,  the  desired  solution  is  actually  self -similar .  In  the 
scaled  variables  £  »  x/t,  n  ■  y/t,  where  t  is  time  and  the  origin  of  ti»e 
coordinate  system  i3  fixed  at  the  wedge  corner,  the  observed  wave  forms  are 
stationary.  It  is  thus  natural  to  seek  a  solution  to  the  ideal  fluid  equa¬ 
tions  in  terms  of  these  similarity  variables. 


Figure 
Fig.  1, 


2.  Triple-point  region  in  ideal  gas  (y  »  1.35)  corresponding  to 
calculated  using  FCT  code  FAST2D. 
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In  the  present  paper  we  treat  the  shock -on-wedge  problem  by  expanding 
the  fluid  variables  and  the  functional  forms  of  the  boundaries  in  power 
series  in  £  and  n.  The  expansion  coefficients  are  found  by  imposing  the 
boundary  conditions  (Rankine-Hugoniot  conditions  on  the  shock,  perfect 
reflection  on  the  wedge  surface). 

Numerous  fluid  dynamics  problems  have  been  solved  by  power  series  tech¬ 
niques,  as  described  by  Van  Dyke.4  Here  we  follow  the  general  approach  he 
advocates:  first  do  the  Icwest-order  terms  "by  hand,"  then  afterward  auto¬ 
mate  the  procedure,  finally  using  sophisticated  mathematical  techniques  to 
find  the  limiting  values  of  critical  numbers  (e.g.,  boundaries  between  two 
reflection  regimes).  Unlike  most  of  the  problems  described  by  Van  Dyke,  the 
algebraic  equations  determining  the  power  series  coefficients  here  are  highly 
nonlinear  and  need  to  be  solved  by  iteration;  the  proper  iteration  scheme  is 
not  obvious,  but  must  be  found  by  experimentation. 

In  the  following  section  we  derive  the  equations  to  be  solved.  Tt\e  next 
section  describes  the  iteration  procedure,  the  program  implementing  it,  and 
the  display.  The  final  section  summarizes  our  conclusions  and  discusses 
extension  of  the  calculation  to  higher  order  and  to  the  more  complicated  Mach 
reflection  cases. 

DERIVATION  OF  EQUATIONS 

In  what  follows  we  use  the  label  "a"  for  quantities  in  the  ambient  gas 
(ahead  of  the  incident  shock);  "b"  behind  the  incident  shock;  "c"  between  the 
contact  surface  and  the  Mach  stem;  and  "d"  for  the  doubly  shocked  region 
between  the  reflected  shock  and  the  contact  surface.  Surfaces  of  gasdynamic 
discontinuities  are  labelled  by  the  two  regions  they  separate,  e.g.  "cd"  for 
the  contact  surface.  The  triple  point  is  labeled  by  the  subscript  ”T, "  and 
the  origin  is  at  the  point  of  attachment  (the  wedge  corner). 
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In  terms  of  the  similarity  variables  £,  n,  the  fluid  equations  become 


3  a  3  o  3u  3v . 

(U-C)  *  (v-n)  _  +  p  (_  +  _)  =  0; 


(1 


P  [(U-C)  ||  +  (v-n)  -f-]  +  -fl  -  0; 


(2 


[(u-5)  IX  +  (v“n)  l£]  +  1?  =  0; 


3n '  3n 


(3 


.  3p  ,  .  3p  ,  3u  3v 

(U-5)  i|  +  (v-n)  +  yp  (-35  +  W 


(4. 


where  Y  is  the  adiabatic  index. 


If  we  define  dimensionless  variables  by 


x  »  £/c,  y  »  n/c  ; 


u  =»  (u-C)c"1  ,  v  =  (v-n)c-1  , 


r  *  p/p  ,  p  =  p/p  f 

d  a 

where 

'2  -  W 


(5) 


(6] 


then  the  fluid  equations  become 


ir"H 


at  n  ■  0  for  all  5.  It  follows  that  we  must  have 


*01  “  RU  *  R  21  *  R31  “  *  *  *  *  0? 


(18) 


U01  3  U11  “  U21  *  U31 


0; 


(19) 


P01  3  P11  3  P  21  3  P  31 


0; 


(20) 


and 


00 


10 


20 


V 


30 


-  0. 


(21) 


To  proceed  further,  we  substitute  Eqs.  (13)— (16)  in  Eqs.  (9)— ( 1 2)  and 
collect  terms  of  like  powers  X1  Y^ ,  equating  their  coefficients  to  zero  order 
by  order.  The  task  of  doing  this  by  conventional  methods  rapidly  becomes 
hopeless,  but  it  can  readily  be  automated.  To  see  what  is  going  on  in  the 
lowest  orders  of  our  "hand"  calculation,  we  expand  using  MACSYMA,  the  sym¬ 
bolic  manipulation  program  developed  at  the  Mathlab  of  MIT.  Hiis  works  up  to 
about  order  N  ■  6,  where  we  define  the  order  of  the  expansion  by  N  * 
max(i+j).  Beyond  that  point  storage  limitations  make  it  necessary  to  use 
"tricks,"  and  eventually  terminate  the  process  entirely. 

when  this  is  done,  certain  patterns  emerge.  We  find  that  we  must  have 


V 


02 


7 1 2  3  7 22 


=  0, 


(22) 


P 03  3  P 1 3  3  P 23 


0. 


(23) 
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Some  of  the  equations  then  vanish  identically.  In  each  of  regions  c  and  d 
there  are  more  variables  than  nontrivial  eauations,  as  shown  by  the  following 
table: 


X'  (U  -  U)  +  Y*  (V  -  V) 


(24) 


Y'  (R  U  -  R  CJ)  *  X'  (R  V  -  R  V); 


P  +  R  (O2  +  V2)  =  P  +  R  (U2  +  V2); 


U2  +  V2  +  Yy  =  U2  +  V2  +  yj  , 


while  on  a  contact  surface 


P  =  P; 


Y'  U  =  X'  V; 


Y'  U  *  X'  V. 


Here  the  variables  in  front  of  and  behind  the  discontinuity  are  denoted 
by  unbarred  and  barred  symbols,  respectively. 


We  know  that  ac  is  nearly  a  straight  line  normal  to  the  wall,  so  the 
representation 


X  -  x0  +  xly  +  x2y *•  +  x 3y 3  +  .  .  .  . 


is  valid.  Substituting  Eq.  (31)  and  the  ambient  conditions  u3  *  v*  *  0  in 
Eq.  (24),  the  condition  for  continuity  of  parallel  velocities  at  ac,  which 
takes  the  form 


X'  (0  -  IT)  +  V  -  V 


0 


(32 


we  find 


Xl  -  Of 


(33) 


1  +  V<01  +  VUX0  +  V21X0  +  *  * 


2  (Xq  +  Uqq  +  U10XQ  +  U20X0  + 


(34) 


etc.  When  the  coefficients  of  Eq.  (31)  ere  known,  substitution  of  Eq.  (31) 


in  Bqs.  (25) -(27)  (with  Y'  ■  1 )  yields  three  conditions  in  each  order  on 

c  c  c  c 

RT .  ,  U.  .  ,  V,  .  ,  P,  .  . 

ID  ij  iD  i  j 


Analogously,  we  write  the  equation  for  the  reflected  shock  bd  as 


Y«  YtX  +  Y2X  +  Y3X  +  .  . 


(35 


Substitution  in  the  condition  for  continuity  of  the  parallel  velocities  at 
bd, 

Ud  -  +  y'  (Vd  -  J°)  -  0,  (36) 


First  we  define  X^,  Y^  by  setting 


*T  -  *1XT  +  *2XT  +  *3XT  * 


XT  *  X0  +  X1*T  +  X2Yt  +  X3Yt  +  *  ‘  * 


Then  we  require 


PC(Xt,Yt)  -  Pd(XT,YT), 


°C<VV  AYt  -  VC(Xt,Yt)  AXt; 


Od(xT,YT)  AYt  -  Vd(XT,YT)  AXt, 


where 


O  -  PC(X,Y)  -  Pd(X,Y)  -  PC(Xt,Yt)  -  Pd(XT,YT) 


*  -  #>T 


determines  the  ratio  Ax_/ Ay_ 


For  our  trial  calculation  we  work  to  order  N  -  5.  Including  the  quanti 


ties  which  vanish  identically,  we  then  have  183  equations  in  183  unknowns, 
which  are  to  be  solved  simultaneously . 

SOLUTION  OF  EQUATIONS 

The  solution  is  obtained  by  iteration.  Several  points  are  important  in 
the  design  of  a  satisfactory  iteration  scheme: 

(i)  Good  values  of  the  quantities  xq,  XT,  YT,  etc.,  can  be 
obtained  from  the  experimental  data  (Fig.  1)  and  used  as  initial 


gu4 


(ii)  The  system  can  be  made  quasilinear  if  we  solve  order  by  order, 
starting  with  N  •  0,  and  using  in  any  given  order  the  previously  obtained 
values  of  all  variables  not  being  solved  for  in  that  order. 

(iii)  To  reduce  the  effect  of  possible  instability  in  parts  of  the 
scheme,  all  variables  are  updated  using  some  form  of  (under) relaxation,  e.g. 
new  value  «  old  value  +  relaxation  factor  x  corrections. 

The  current  version  of  the  code  works  as  follows: 

(i)  For  N  »  0,  use  the  X-  and  Y-  independent  terms  of  Eqs.  (9)-(12), 


together  with  the  zeroth-order  form  of  Eqs.  (25) -(27)  and  Bq.  (43)  to  update 


c  c  c  c  c  c  c 
R  ,  R  ,  U  ,  U  ,  V  ,  P  ,  P 


00  10  00  10  01  00 
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(ii)  Do  the  same  thing  [using  Eq.  (44)  instead  of  (43)]  for  R  , 


Pd  ; 
10 


00 


(iii)  Solve  Bqs.  (9)-(12)  and  Eqs.  (25)-(27)  on  ac  order  by  order  for 
N  »  1  to  5  to  obtain  the  region-c  variables; 

(iv)  Solve  Eqs.  (9) -(12)  and  Eqs.  (24) -(27)  on  bd  order  by  order  for 
N  ■  1  to  5  to  obtain  the  region-d  variables  plus  YN+1 ; 

(v)  Iterate  Eqs.  (40)-(42)  together  with  the  successive  orders  of 
Bq.  (32)  [e.g.,  Bq.  (34)  for  X2)  to  solve  for  X^,  YT,  and  X^  i  ■  0,  1, 
....,  6.  (To  a  good  approximation,  X,  *  0  for  all  i  *  0,  2.) 


This  scheme  is  repeated  until  no  further  change  (to  some  preset  tolerance)  in 
the  variables  occurs*  The  program,  written  in  Fortran,  runs  on  a  VAX  11/780 
at  about  one  second  per  iteration. 

Results  are  most  conveniently  displayed  as  contour  plots  in  pressure  (or 
density) .  Although  it  is  possible  to  drive  the  plotting  pen  directly  using 
the  exact  formula  P(X,Y),  it  is  easier  to  declare  a  very  large  array  (e.g., 
500*500),  fill  it  with  pressures  calculated  at  every  X,  Y,  and  use  a  standard 
contour  plotting  package  on  this. 

CONCLUSION 

The  program  describe'  above  is  still  under  development,  and  no  results 
have  yet  been  generated.  For  this  sample  problem,  it  seems  straightforward 
to  carry  the  method  to  a  successful  solution.  We  have  encountered,  and 
apparently  overcome,  two  types  of  difficulties:  determining  the  formulation 
for  the  problem,  and  solving  the  resulting  set  of  equations.  Only  if  both 
stages  are  handled  correctly  will  useful  answers  results.  There  remains  the 
(programming)  task  of  automating  the  solution,  so  as  to  work  to  arbitrarily 
high  order  N.  This  is  of  interest  chiefly  in  connection  with  studying  the 
behavior  of  the  roll-up  in  the  contact  surface. 

Finally,  extension  to  other  forms  of  Mach  reflection  is  of  interest. 
Attached  double  (or  triple)  Mach  reflection  presents  no  problem  in  principle, 
and  can  be  handled  using  the  techniques  described  here.  Complex  Mach 
reflection  and  detached  shocks  are  less  clear.  At  present  we  do  not  know  how 
to  formulate  the  problem  in  either  situation  so  as  to  make  it  well-posed. 

This  will  be  addressed  in  future  work. 
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Olcy  Actn  Mail  Ser  Sec  L  Vortman 

Olcy  Attn  Mail  Ser  Sec  A  Chaban 

Olcy  Attn  Mail  Ser  Sec  L  Hill 

Science  Applications,  Inc. 

P  0  Box  2351 
La  Jolla,  CA  92038 

Olcy  Attn  Technical  Library 

Science  Applications,  Inc. 

1257  Tasman  Drive 
Sunnyvale,  CA  94086 
Olcy  Attn  J  Dishon 

Science  Applications,  Inc. 

2450  Washington  Avenue 
San  Leabdro,  CA  94577 
Olcy  Attn  2  Maxwell 

Olcy  Attn  D  Bernstein 

Science  Applications ,  Inc. 

?  0  Box  1303 
McLean,  VA  22102 

Olcy  Attn  M  Knasel 

Olcy  Attn  B  Chambers  III 

Olcy  Attn  R  Sievers 

Olcy  Actn  J  Cockayne 

Southwest  Research  Institute 
P  0  Drawer  28510 
San  Antonio,  TX  78284 
Olcy  Attn  W  Baker 
Olcy  Attn  A  Wenzel 

SRI  International 
333  Ravensvood  Avenue 
Menlo  Park,  CA  94025 

Olcy  Actn  G  Abra'namson 

Systems,  Science  &  Software,  Inc. 

?  0  Box  1620 
La  Jolla,  CA  92038 
Olcy  Attn  Library 
Olcy  Attn  D  Crine 
Olcy  Art r.  T  Riney 
Olcy  Attn  F  Pvatt 
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Teledyne  Brown  Engineering 
Cummings  Research  Park 
Huntsville,  AL  35807 

Olcy  Attn  J  Ravenscraft 

Terra  Tek,  Inc. 

420  Wakara  Way 
Salt  Lake  City,  UT  84108 
Olcy  Attn  Library 
Olcy  Attn  S  Green 
Olcy  Attn  A  Jones 

Tetra  Tech,  Inc. 

630  N  Rosemead  Blvd. 

Pasadena,  CA  91107 
Olcy  Attn  L  Hwang 
Olcy  Attn  Library 

TRW  Defense  &  Space  Sys  Group 
One  Space  Park 
Redondo  Beach,  CA  90278 
Olcy  Attn  I  Alber 
Olcy  Attn  Tech  Infor  Ctr 
02cy  Attn  N  Lipner 
Olcy  Attn  P  Bhutta 
Olcy  Attn  D  Baer 
Olcy  Attn  R  Plebuch 

TRW  Defense  &  Space  Sys  Group 
P  0  Box  1310 

Sar.  Bernardino,  CA  92402 
Olcy  Attn  E  Wcng 
Olcy  Attn  P  Dai 

Universal  Analytics,  Inc. 

7740  W  Manchester  Blvd 
Playa  Del  Rey,  CA  90291 
Olcy  Attn  E  Field 

Veidlinger  Assoc.,  Consulting  Eng 
110  E  59th  Street 
New  York,  NY  10022 
Olcy  Attn  M  Baron 

Weidlinger  Assoc.,  Consulting  Eng 
3000  Sand  Hill  Road 
Menlo  Park,  CA  94025 

Olcy  Attn  J  Isenberg 

Westinghouse  Electric  Corp. 

Marine  Division 
Her.dy  Avenue 
Sunnyvale,  CA  94088 
Olcy  Attn  W  Volz 
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